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ABSTRACT 

_ ^ 

'  A  review  is  presented  of  remote  sensing  methods  and  associated  inverse  procedures 
for  obtaining  the  acoustic  properties  of  sediments  in  shallow  water.  The  properties  of 
interest  are  the  vertical  profiles  of  density,  compressional  wave  velocity  and  attenuation, 
and  shear  wave  velocity  and  attenuation.  For  compressional  waves,  the  remotely  sensed 
information  required  is  the  modal  reflection  coefficient  of  the  seabed,  usually  expressed 
in  terms  of  the  depth-dependent  Green’s  function  at  specific  locations  in  the  water 
column.  For  shear  waves,  information  about  the  Scholte  and  Rayleigh  interface  waves 
on  the  seabed  is  necessary.  Also  reviewed  are  methods  for  remote  sensing  of  the  Biot 
geophysical  parameters,  shear  modulus,  permeability,  and  porosity,  from  which 
geoacoustic  properties  can  be  derived. 
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I.  INTRODUCTION 

In  a  previous  report,1  the  geoacoustic  and  geophysical  properties  of  sediments  were 
discussed  in  their  relation  to  shallow  water  acoustics  at  low  frequencies  (20  to  500  kHz). 
The  purpose  of  the  present  report  is  to  review  current  methods  of  remote  sensing  and 
inversion  of  measurement  data  in  shallow  water  to  obtain  these  properties.  We  shall 
consider  the  acoustic  spectrum  in  three  different  frequency  regions:  (a)  £  100  Hz,  (b) 
10-100  Hz,  and  (c)  1  to  10  Hz. 

At  frequencies  2: 100  Hz,  the  reflectivity  and  scattering  behavior  of  the  seabed  are 
used  to  determine  its  geoacoustic  properties  where  the  underlying  sediment  is  assumed  to 
be  a  viscoelastic  fluid2  rather  than  a  simply  lossless  fluid.  Hydrophones,  spaced 
vertically  and  in  range  at  planned  locations,  allow  the  processing  of  measurement  data 
from  either  monochromatic  cw  or  impulse  sources  to  yield  the  required  profiles.  In  this 
approximation,  generation  of  shear  waves  is  treated  as  an  added  loss  to  the  compressional 
wave  attenuation. 

An  important  problem  in  shallow  water  acoustics  is  the  determination  of  those 
geophysical  properties  of  the  seabed,  sedimentary  layers,  and  acoustic  basement  that 
account  for  the  behavior  of  acoustic  waves  at  frequencies  below  about  50  Hz.  Although 
at  these  frequencies  acoustic  energy  can  penetrate  deeply,  there  is  usually  little  or  no 
direct  access  to  the  geological  structure  for  determining  the  actual  geoacoustic  properties. 
Remote  sensing  techniques  must  therefore  be  developed  from  which  density  profiles  and 
compressional  and  shear  wave  profiles  of  velocity  and  attenuation  can  be  inferred.  Using 
this  information,  the  acoustic  wave  equation  can  be  solved. 

In  the  frequency  range  below  about  50  Hz,  the  poroviscous  properties  of  the 
sediments  must  be  treated  using  Biot  theory.  The  consequences  of  Biot  theory  in  this 
case  have  been  studied  by  Stoll3  and  Yamamoto,4  who  have  also  applied  inverse  theory 
and  methods  to  the  problem.  The  two  researchers  are  particularly  interested  in 
permeability,  porosity,  viscosity,  and  density.  These  quantities  determine  the 
poroviscous  frequency  number  and  the  relaxation  frequency  below  which  the  attenuation 
falls  off  approximately  as  the  square  of  the  frequency. 

Below  the  acoustic  cutoff  of  the  shallow  water  sound  channel,  at  approximately  1  to 
10  Hz,  energy  is  transmitted  by  interface  waves  along  the  bottom.5  Knowledge  of  Scholte 
and  Rayle.'gh  wave  traces  on  the  seabed  allows  the  determination  of  the  shear  wave 
properties  of  the  sedimentary  and  basement  layers.  Some  previously  developed  seismic 
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techniques  can  also  be  employed  for  this  purpose.  The  properties  of  the  longer  Rayleigh 
waves  are  used  to  determine  the  thickness  of  the  sediment  layer  as  well  as  the  material 
properties  of  the  acoustic  basement.  Furthermore,  the  dispersion  of  interface  waves  and 
the  modal  excitation  of  boundary  and  related  shear  waves,  along  with  time  signals 
stacked  in  range  (synthetic  seismograms),  allow  the  determination  of  the  seismoacoustic 
parameters  of  the  sediment  layer.  These  determinations  are  usually  carried  out  by  fast 
forward  methods. 

A  forward  method  starts  with  a  model  of  the  scattering  and  transmission  physics, 
usually  in  the  form  of  acoustic  pressure  measurements  and  propagation  loss  runs  versus 
range  and  depth.  As  model  parameter  values  are  varied,  computer  solutions  are  iterated 
until  a  best  fit  of  model  to  data  is  obtained  in  some  defined  sense.  Uniqueness  and 
consistency  of  the  final  solution  are  not  guaranteed.  To  rectify  this,  a  theory  of  direct 
inversion  has  been  developed  that  allows  the  treatment  of  uniqueness  and  optimal 
experimental  design  in  the  presence  of  interfering  noise.  Theoretically,  with  this 
technique  one  can  use  as  a  starting  point  the  same  model  as  in  the  forward  case,  though 
the  experimental  Technique  might  require  a  different  space  and  time  configuration  of 
measurement  data. 

The  problem  then  becomes  one  of  designing  efficient  experiments  using  sources 
whose  properties  are  known  with  respect  to  frequency  content,  energy  content,  and  the 
time  domain  signatures.  If  appropriate  receivers  are  placed  properly  in  range  and  depth, 
including  bottom  mounting,  the  data  can  be  inverted  to  obtain  the  acoustic  description  of 
the  sedimentary  stratification  and  the  underlying  basement  with  some  guarantee  of 
uniqueness. 

The  problem  of  inverse  experimental  design  and  data  analysis  as  applied  to 
scattering  experiments  has  had  a  long  hL  try  in  the  solution  of  boundary  value  problems 
associated  with  the  wave  equations  of  quantum  mechanics,  nuclear  physics, 
electromagnetics,  seismology,  and  acoustics.  The  Appendix  is  a  summary  of  inverse 
theory  for  the  linear  and  nonlinear  cases  of  shear  modulus  profiling  adapted  from 
Yamamoto  and  Torii.6  In  this  report,  I  review  remote  sensing  methods  and  inversions 
presently  used  to  find  the  acoustic  properties  of  the  bottom  and  subbottom  as  they  affect 
shallow  water  propagation.  However,  I  will  consider  only  those  cases  where  all  three 
quantities  -  density,  velocity,  and  attenuation  -  are  determined  or  where  geophysical 
parameters  are  measured  that  allow  their  determination.  These  geophysical  parameters 
include  permeability,  the  complex  shear  modulus,  and  porosity. 
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a  COMPRESSIONAL  WAVE  PROPERTIES 
A.  Introduction 

Detailed  reviews  exist  for  the  problem  of  remote  sensing  and  inverse  methods  for 
compressional  wave  properties  of  shallow  ocean  sediments.7'11  Here  I  outline  only  the 

A 

major  points. 

The  basic  remote  sensing  information  required  is  the  aperture-  and  band-limited 
values  of  reflected  acoustic  pressure  waves  as  a  function  of  depth  in  the  water  column.  In 
deep  water,  the  theoretical  effort  centers  on  a  model  of  plane  waves  reflected  by  a  one¬ 
dimensional  inhomogeneous  fluid  layer.12  In  shallow  water,  the  reflected  plane  wave  is 
treated  in  modal  form  as  the  depth-dependent  Green’s  function  versus  horizontal  wave 
number.13 

Compressional  wave  modeling  in  shallow  water  sediments  usually  assumes  that  the 
sediment  is  a  horizontally  layered  viscoelastic  fluid.  Shear  properties  and  Biot  theory 
refinements  are  thus  ignored  for  a  first  approximation.  On  this  basis,  both  forward  and 
direct  inverse  methods  have  been  used  to  model  the  density  and  sound  velocity  profiles 
with  depth.  If  the  sediment  attenuation  coefficient  is  linear  with  frequency,  then  the 
attenuation  profile  can  also  be  derived  for  this  fluid  medium  model. 

Forward  methods,  such  as  inferring  bottom  loss  from  propagation  loss  versus  range 
experiments,  have  the  disadvantage  of  possible  incompleteness  of  model  parameters. 
Direct  inverse  methods  can  test  the  uniqueness  of  models  and  ern  be  used  efficiently  for 
designing  appropriate  experiments. 


B.  Inverse  Methods 

1 .  Plane  Wave  Reflection  Coefficient  versus  Angle  of  Incidence 

The  plane  wave  reflection  coefficient  of  the  ocean  bottom  shows  resonances  both  as 
a  function  of  frequency  and  angle  of  jcidence,  depending  on  the  properties  of  the  ocean 
bottom.  This  is  valid  for  any  water  depth  since  its  input  is  the  plane  wave  reflection 
coefficient  of  the  bottom,  which  is  independent  of  the  water  depth.7,8,11,12,14,15  However, 
there  are  practical  limitations  in  shallow  water. 


*!  have  seen  only  one  reasonably  complete  treatment  of  the  actual  inversion  itself  that  start¬ 
ed  with  at-sea  measurements  for  compressional  wave  properties.7 
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Catholic  University  (CU)  method  of  inversion 

The  CU  method7,11,15  analyzes  how  the  bottom  reflection  coefficient  is  affected  by 
layer  resonances  that  are  carried  by  acoustic  signals  reflected  from  the  ocean  floor.  The 
method  exploits  the  fact  that  measured  resonance  parameters  (locations,  heights,  and 
widths)  are  related  to  the  layered-medium  parameters  (densities,  layer  thicknesses,  and 
sound  velocities  of  the  layers)  and  hence  provide  a  solution  to  the  inverse  problem.  The 
model  assumes  a  three-layer  system  -  water,  sediment,  and  basement.  See  Figure  1, 
where  d  is  the  sediment  layer  thickness,  p  is  the  density,  c  is  the  sound  velocity,  and  0  is 
the  angle  of  incidence. 


Figure  1.  Fluid  layer  embedded  between  two  different  fluids  with  incident  and 
transmitted  sound  wave  (from  Ref.  7). 

If  the  square  of  the  reflection  coefficient  is  plotted  on  a  three-dimensional  graph 
versus  the  frequency  (/)  layer  thickness  product  fd  and  incident  angle  03  (from  water  to 
sediment  angle),  one  gets  a  surface  with  a  series  of  ridges  and  valleys.  By  plotting  this 
upside  down  as  the  square  of  the  transmission  coefficient,  we  emphasize  the  sharp  ridges 
descriptive  of  resonance-type  processes.  Nagl  ct  al.7, 11,15  proceed  by  interpreting  the 
scattering  process  as  a  Breit-Wigner  process  of  quantum  physics  for  which  the  inverse 
solution  already  exists. 

The  experimentally  measured  resonance  parameters  are  then  related  to  the  acoustic 
properties  of  the  reflecting  medium.  By  measuring  the  critical  angle,  the  minimum  value 
of  the  reflection  coefficient,  the  spacing  and  width  of  the  fd  resonances  at  some  angle, 
and  the  position  and  width  of  an  angle  resonance  at  a  particular  frequency,  all  the 
parameters  describing  the  reflecting  medium  for  a  simple  case  —  c  j,  c2,  pj,  p2,  d  —  can 
be  deduced.  Because  the  positions  and  widths  of  the  angular  resonances  vary  with  angle, 
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it  should  be  possible  to  determine  a  much  larger  set  of  medium  parameters  than  for  the 
simple  case  considered  here.  However,  the  method  has  not  been  exploited  for  either  deep 
or  shallow  water. 

D.J.  Thomson14  has  also  investigated  the  problem  of  determining  the  material 
properties  of  the  seabed,  i  e.,  the  density  and  sound  speed  profiles  as  a  function  of  depth 
from  an  aperture*  and  band-limited  knowledge  of  the  acoustic  reflection  coefficient.  He 
uses  a  direct  inverse  solution  method  from  Candel  et  al.9  This  method  is  based  on  the 
scattering  of  plane  waves  by  a  one-dimensional  inhomogeneous  medium.  Thomson’s 
use  of  the  Bom  approximation  provides  him  a  high-frequency  approximation,  whereas 
the  method  of  Candel  et  al.  uses  the  low-frequency  information  to  estimate  the  gradients 
for  sound  speed  and  density.  Thomson  presents  an  equation  for  the  reflection  coefficient 
that  supports  a  noniterative  method  to  deduce  the  admittance  versus  depth  profile  of  an 
inhomogeneous  medium  from  the  impulse  response.  This  equation  is  a  nonlinear  Fourier 
transform  of  the  logarithmic  derivative  of  the  local  admittance.  Inversion  of  the  integral 
transform  enables  the  recovery  of  the  admittance  profile  using  the  numerical  integration 
of  two  first-order  differential  equations  that  require  a  single  impulse  response  of  the 
medium  as  reflection  data.  He  concludes  that  reconstruction  of  both  the  density  profile  p 
(z)  and  the  sound  velocity  profile  c(z)  should  be  no  more  difficult  than  determining  one  of 
these  quantities.  It  is  necessary  to  have  reflection  responses  for  two  different  probing 
directions,  and  four  first-order  differential  equations  have  to  be  integrated  instead  of  two. 

Finally,  Thomson  points  out  that  the  derivation  of  the  solution  requires  that  the 
derivative  of  the  logarithm  of  the  vertical  admittance  in  the  vertical  direction  be 
independent  of  frequency.  This  requirement  is  satisfied  for  nonabsorbing  media.  It  is 
also  satisfied  for  absorbing  media  if  the  absorption  depends  linearly  on  frequency.  There 
is  evidence  that  this  is  not  true  for  sediment  absorption  at  the  low  frequencies  below  the 
relaxation  frequency  for  sands  and  silts  (see  Section  IV).  Furthermore,  the  method  does 
not  apply  to  geoacoustic  models  that  support  multiple  subbottom  reflections. 

2.  Depth-Dependent  Green’s  Function 

The  reflection  coefficient  can,  in  principle,  be  obtained  from  the  depth-dependent 
Green’s  function.  A  method  for  inverting  acoustic  pressure  measurements  taken  in 
shallow  water  using  the  depth-dependent  Green’s  function  has  been  developed  by  Rajan 
and  his  associates  at  the  Woods  Hole  Oceanographic  Institution  (WHOI).7,10,13,16-18  The 
required  input  data  are  the  magnitude  and  phase  of  the  pressure  field  versus  range  with  a 
monochromatic  (cw)  point  source. 
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Synthetic  aperture  array  data  are  obtained  by  towing  a  cw  source  away  from  moored 
receivers.  The  data  are  Hankel-transformed  numerically  to  get  the  depth-dependent 
Green’s  function  versus  horizontal  wave  number.  In  terms  of  normal  mode  theory,  the 
Green’s  function  contains  information  about  the  nature  of  the  discrete  and  continuous 
modal  spectra  and  the  plane-wave  reflection  coefficients  of  the  waveguide  boundaries. 
Specificall*’,  the  Green’s  function  contains  sharp  peaks  at  horizontal  wave  numbers 
corresponding  to  the  eigenvalues  of  trapped  and  virtual  modes  excited  in  the  waveguide. 
The  positions  and  magnitudes  of  these  modal  peaks  can  be  sensitive  to  the  acoustic 
properties  of  the  bottom. 

A  geoacoustic  model  is  then  determined  by  matching  the  bottom  model  prediction 
of  the  depth-dependent  Green’s  function  to  the  experimental  Green’s  function, 
particularly  in  the  positions  of  the  modal  peaks.  Forward  model  iterative  techniques  and 
perturbative  inverse  techniques  have  both  been  applied  successfully  to  this  inverse 
problem. 

The  WHOI  method  involves  an  exercise  in  which  a  cw  source  is  towed  at  a  fixed 
depth  away  from  two  moored  receivers  over  an  aperture  that  extends  from  zero  range  to 
several  kilometers.  The  configuration  is  shown  in  Figure  2,  where  Z  is  depth.  The 
receivers  quadrature-demodulate  the  signal  by  beating  the  signal  down  to  0  Hz  to  remove 
the  harmonic  time  dependence  and  then  record  the  real  and  imaginary  aspects  of  the 
spatial  part  of  the  acoustic  pressure  field  p(r). 

The  forward  model  iteration  inversion  consists  of  computing  the  theoretical  Green’s 
function  for  various  values  of  the  bottom  parameters  and  finding  in  the  positions  and 
relative  amplitudes  of  the  modal  peaks  the  parameter  set  that  pro  /ides  the  best  agreement 
with  the  experimental  Green’s  function.  Comparisons  are  also  made  between  the 
measured  and  the  theoretically  computed  pressure  field  magnitudes,  so  that  a  consistent, 
balanced  fit  is  obtained  in  both  domains. 

For  the  perturbative  inverse  technique,13,17,18  the  value  of  the  location  of  the 
trapped-mode  eigenvalues  is  the  data  used  to  invert  for  the  difference  Ac(z)  between  the 
actual  sound  speed  profile  of  the  bottom  c(z)  and  that  of  a  plausible  reference  model. 
Specifically,  one  locates  the  mode  peaks  and  subtracts  the  mode  peak  locations  to  form 
Akn  for  each  of  the  modes  and  then  uses  these  as  input  to  an  integral  equation  for  Ac  (z). 
Tne  symbol  kn  stands  for  the  eigenvalues  of  the  horizontal  wavenumber  k , 

The  equation  for  the  Green’s  function  is  obtained  as  a  straightforward  application  of 
first-order  perturbation  theory.  It  is  readily  extended  to  find  attenuation  and  density 
profiles.13,17,18 
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A  demonstration  of  the  technique  starts  with  Figure;-  3  and  4,  which  show  the 
pressure  fields  measured  for  a  sandy  bottom  over  a  range-aperture  of  1325  m  in 
Nantucket  Sound,  Massachusetts,  May  1984. 
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Figure  2.  Typical  shallow  water  experimental  setup  (from  Ref.  7). 


Figure  3.  Comparison  of  theoretical  and  experimental  pressure  field  magnitudes  and 
phases  at  140  Hz,  using  forward  model  iteration  (from  Ref.  7). 
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Figure  4.  Comparison  of  theoretical  and  experimental  pressure  field  magnitudes  and 
phases  at  220  Hz,  using  forward  model  iteration  (from  Ref.  7). 

The  magnitudes  show  strong  spatial  interference  patterns.  These  patterns  arise  from 
the  coherent  combination  of  two  normal  modes  of  comparable  strength  at  220  Hz  and 
one  dominant  and  one  weak  mode  at  140  Hz.  The  decay  with  range  of  the  weaker  mode 
at  140  Hz  is  also  evident  in  the  data.  In  addition,  an  increase  of  v'ater  depth  from  13.9  to 
14.6  m  at  a  range  of  about  600  m  has  a  pronounced  effect  on  the  220  Hz  modal  cycle 
distance,  which  changes  in  that  vicinity  from  about  125  to  150  m.  The  results  of  the 
forward  model  iteration  inversion  are  shown  in  the  theoretical  curves  for  p(r)  in 
Figures  3  and  4.  Excellent  agreement  is  obtained  for  the  220  Hz  data,  but  only  fair 
agreement  for  the  140  Hz  data.  The  likely  reason  for  such  a  failure  is  probably  the 
inherent  simplicity  of  the  bottom  model,  where  only  one  or  two  parameters  were  varied 
at  a  time.  A  more  general  approach  such  as  direct  inverse  theory  is  required.  Here  the 
bottom  is  not  constrained  to  a  few  isovelocity  or  simple  gradient  layers. 

In  Figure  5,  the  pressure  field  generated  by  using  the  perturbative  inverse  scheme  is 
shown  for  the  140  Hz  data. 

The  agreement  between  theory  and  measured  data  is  also  excellent  at  220  Hz,  which 
is  not  shown.  Note  that  this  agreement  was  obtained  after  only  three  iterations  of  an 
inversion  program  that  used  a  crude  Pekeris  model  as  a  background.  This  demonstrates 
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the  efficiency  and  robustness  of  the  direct  inverse  approach  as  compared  with  standard 
forward  modeling. 


Figure  5.  Comparison  of  theoretical  and  experimental  pressure  field  magnitudes  at 
140  Hz,  using  perturbative  inverse  technique  (from  Ref.  7). 


Figure  6.  Geoacoustic 
models  of  the  bottom 
produced  by  forward 
modeling  and  by  pertur¬ 
bative  inverse  technique 
(from  Ref.  7). 


Figure  6  shows  the  geoacoustic  models  produced  by  the  two  techniques.  The  high 
near-surface  sound  velocity  gradients  produced  by  both  models  are  consistent  with 
Hamilton’s  results  for  sandy  bottoms.19  The  perturbative  inverse  technique  also 
produces  a  continuously  variable  gradient  profile  until  constrained  at  a  sediment  depth  of 
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10  m.  This  allows  the  fit  of  data  that  cannot  be  handled  by  simpler  models  such  as  the 
constant  gradient  layer  used  in  the  forward  approach.  The  symbol  a  in  Figure  6  stands 
for  attenuation. 

Rajan  reports  that  the  low-frequency  inverse  problem  becomes  linear  as  frequency 
is  decreased  and  that  linearization  of  high-frequency  inverse  problems  can  be  obtained 
from  this  solution. 

3.  NORDA  Method 

The  NORDA  method7,21-23  involves  two  basic  concepts  in  terms  of  the  depth- 
dependent  Green’s  function.  First,  the  magnitude  of  the  beam  response  of  a  horizontal 
array  is  approximately  related  to  the  magnitude  of  the  depth-dependent  Green’s  function. 
Specifically,  the  beam  response  will  exhibit  maxima  at  angles  corresponding  to  the 
propagation  angles  of  trapped  and  virtual  modes.  The  angular  resolution  is  in  principle 
inversely  proportional  to  the  range  aperture.  In  practice,  it  is  still  good  for  short  array 
lengths  and  small  numbers  of  receivers. 

Second,  the  method  uses  the  idea  that  for  lossy  bottoms  the  modal  attenuation  with 
range  increases  with  mode  number.  This  assumes  that  for  bottom  structures  consisting  of 
isovelocity  layers  with  monotonically  increasing  velocity,  the  higher  modes,  with  larger 
grazing  angles,  penetrate  into  the  bottom  more  deeply  and  suffer  greater  attenuation.  As 
a  result,  the  source-receiver  array  offset  is  increased  and  the  beam  response  is  dominated 
by  lower  modes  that  interact  with  shallower  portions  of  the  bottom.  Since  a  given  mode 
with  eigenangle  0m  can  be  viewed  as  being  trapped  by  the  reflecting  horizon  with  lowest 
sound  velocity  cm  such  that  cm  £  d cos  0m,  the  variable-offset  beamformer  response  can 
be  used  to  infer  sound  velocities  in  the  bottom. 

As  the  source-receiver  offset  is  progressively  increased,  beamformed  responses 
from  deeply  penetrating  modes  also  progressively  disappear  owing  to  high  in-situ 
absorption.  Conversely,  beamformed  peaks  corresponding  to  progressively  shallower 
layers  dominate  the  response  in  turn.  Previous  work23  has  shown  that  it  is  usually 
possible  to  identify  compressional  and  shear  wave  velocity  in  a  range  of  multilayered 
bottoms  by  moving  the  array  farther  away  from  the  source  and  systematically  isolating 
the  critical  angles  of  each  layer. 
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HI.  SHEAR  WAVE  PROPERTIES* 

A.  Introduction 

According  to  Akal  and  Jensen,24  the  excitation  of  shear  waves  in  the  bottom  often 
becomes  the  dominant  loss  mechanism  for  waterborne  sound  in  shallow  water.  These 
authors  claim  that  a  realistic  model  of  the  ocean  bottom  is  a  viscoelastic  solid  described 
by  compressional  and  shear  wave  velocities,  attenuations,  and  densities.  Density  and 
compressional  wave  speed  can  be  found  by  direct  methods,  but  shear  wave  properties  are 
difficult  to  measure:  their  attenuation  is  high,  and  it  is  difficult  to  generate  waves  with 
predominantly  transverse  particle  motions.  However,  the  required  parameters  can  be 
determined  from  measurements  of  interface  waves  on  the  ocean  bottom.  While  shear 
properties  of  core  samples  can  be  measured  in  the  lab,  indirect  determinations  rely  on  the 
measurement  of  the  boundary  wave  parameters  of  sediment-water,  sediment-basement, 
and  basement-water  surfaces.  There  are  three  kinds  of  closely  related  interface  waves: 
Rayleigh,  Stoneley,  and  Scholte.  Rayleigh  waves  propagate  along  the  surface  of  a  free 
solid.  Stoneley  waves  are  bound  to  the  interface  between  two  solids,  i.e.,  sediment  and 
acoustic  basement.  Scholte  waves  propagate  along  a  liquid-solid  interface. 

Rauch5  has  reviewed  the  role  of  interface  waves  in  ocean  seismoacoustics.  He 
discusses  the  importance  of  the  low  velocity  Scholte  waves  for  seismic  sensing  of 
infrasound  in  coastal  waters.  Such  waves  are  below  waveguide  cutoff.  In  seismics,  the 
frequency  0.2  Hz  separates  short  period  tremors  from  medium  and  long  period  earth 
movements.  For  acoustic  studies,  30  Hz  divides  infrasonics  from  sonics.  In  the 
following  sections  I  concentrate  on  the  problem  of  deducing  the  shear  properties  of 
sediments  from  interface  waves. 

B.  Seismic  Interface  Waves 

The  usual  method  of  exciting  seismic  interface  waves  in  shallow  water  is  to  use 
broadband  sources  such  as  charges,  airguns,  water  guns,  and  sledges.  There  are  cw 
vibration  exciters  for  individual  frequencies  and  also  for  FM  pulses.  The  usual  receivers 
for  interface  waves  are  three-axis  geophones  mounted  on  the  bottom  along  with  a  nearby 
hydrophone.  This  combination  is  called  an  ocean  bottom  seismometer  (OBS).  Because 

)|( 

The  figures  in  this  chapter  have  been  taken  from  different  sources  where  symbols  may  have 
different,  but  clear,  meanings. 
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the  longer  wavelengths  of  lower  frequencies  penetrate  deeper,  the  depth-dependence  of 
shear  wave  properties  can  be  found  by  varying  frequency.  The  low-frequency  content  of 
broadband  signals  is  also  used  for  this  purpose. 

The  three  types  of  interface  waves  have  the  following  properties  in  common:  (1) 
they  are  damped  exponentially  with  distance  from  the  interface  in  both  directions,  (2) 
their  particle  motion  is  elliptically  polarized  in  the  depth-iange  plane,  (3)  their  velocity  is 
closely  related  to  the  shear  velocity  of  the  solids  involved,  and  (4)  there  is  no  low- 
frequency  cutoff.  Many  acoustic  measurements  in  shallow  water  have  shown  that 
infrasound  is  transmitted  below  the  acoustic  cutoff  frequency  of  the  water  duct  as  Scholte 
interface  waves. 

Rauch  and  Schmalfeldt25  have  measured  the  characteristics  of  interface  waves 
generated  by  0.18-kg  TNT  shots  at  a  range  of  Ik  -  in  a  water  depth  of  25 m.  From  an 
OBS  they  plotted  sensor  output  versus  time,  incluoing  particle  velocities  from  a  three- 
axis  geophone  and  the  pressure  output  from  the  associated  hydrophone  (see  plots  in 
Figure  7,  which  also  shows  particle  velocity  hodographs).  A  pronounced  long-lasting 
and  low-frequency  interface  wavelet  is  seen  to  arrive  much  later  than  the  high-frequency 
onset  of  the  water  wave.  The  recorded  signals  and  hodographs  show  that  the  interface 
wave  is  detected  only  by  the  vertical  and  radial  geophones.  There  is  little  or  no  energy  in 
the  transverse  direction. 

C.  Dispersion  of  Interface  Waves 

Instead  of  using  shots,  Holt  et  al.26  used  a  bottom-mounted  vibration  exciter  as  a  cw 
source  to  generate  interface  waves  at  individual  selected  frequencies.  They  then 
determined  the  dispersion  of  Scholte-  and  Rayleigh-type  waves  and  thereby  deduced  the 
shear  modulus  profile  of  the  seafloor. 

They  studied  the  propagation  characteristics  of  the  cw  pulses  in  detail.  The 
received  pulses  had  sound  velocities  in  the  range  100-200  m/s  that  are  characteristic  of 
Scholte  waves.  While  Scholte  waves  could  also  be  studied  using  hydrophones  near  the 
bottom,  these  are  less  sensitive  to  signals  and  more  sensitive  to  ambient  noise  than 
geophones. 
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Figure  7.  (a)  A  plot  of  shot-generated  time  series  from  three  orthogonal  geophones 
and  one  hydrophone,  (b)  Hodographs  of  the  particle  velocity  for  a  2-s  time 
window  (from  Ref.  25).  Details  may  be  obtained  from  the  authors.  For  our 
purposes,  these  should  be  understood  as  scaled  amplitudes. 
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The  authors  first  measured  the  phase  velocity  as  a  function  of  frequency  between 
5  Hz  (limited  oy  the  vibrator)  and  100  Hz  (limited  by  the  seismometers).  They  then 
found  a  best  fit  to  the  data  as  shown  in  Figure  8.  With  this  best  fit,  they  then  calculated 
the  group  velocity  using  the  standard  formula 


Vt(f)~ 


VpC f) 

f  dVp(f)  • 
VpC f)  df 


(1) 


where  /  is  the  frequency  in  hertz,  Vp  is  the  phase  velocity  in  meters  per  second,  and  Vg 
is  the  group  velocity  in  meters  per  second. 


Frequency  [Hz] 

Figure  8.  Measured  phase  velocity  from  continuous  wave  method  (■)  and  filtered  pulse 
technique.  (Hatched  areas :  width  given  by  the  bandwidth  of  the  filter; 
height  given  by  the  phase  variations  across  the  pulse.  Average  values  are 
denoted  by  horizontal  lines.)  (•)  denotes  measured  phase  velocity  at  certain 
given  frequencies  in  the  pulses.  The  group  velocity  ( broken  line)  is  calcu¬ 
lated  from  the  fully  drawn  phase  velocity  curve  (best  fit)  (from  Ref.  26). 
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The  inverse  problem  is  to  find  the  shear  velocity  of  the  sediment  layer  as  a  function 
of  depth  using  the  geoacoustic  model  of  Figure  9.  The  water  column  and  solid  substrate 
are  homogeneous  and  semi-infinite,  but  the  density  and  elastic  parameters  of  the 
sediment  can  vary  arbitrarily  with  depth.  The  depth-separated  wave  equation  of  the 
inhomogeneous  sediment  layer  is  solved  by  numerical  integration.27 


Figure  9.  The  physical  model  (from  Ref.  26). 

Holt  et  al.  solved  the  wave  system  for  the  model  shown  in  Figure  9  as  follows.  If 
A0  is  the  upgoing  compressional  wave  in  the  water  and  there  are  two  downgoing  waves 
in  the  solid  substrate  (a  compressional  wave)  and  B2  (a  shear  wave),  then  setting 
boundary  conditions  at  the  two  sediment  interfaces  allows  the  determination  of  the 
dispersion  equation.28  The  boundary  conditions  require  that  the  stress-displacement 
vector  is  continuous  at  the  interfaces.  In  Figure  9,  p  is  the  density,  v,  is  the  shear  wave 
velocity,  v_,  is  the  compressional  wave  velocity,  c0  is  the  sound  speed  in  water,  and  z  is 
depth. 

Using  this  theoretical  model,  the  sediment  parameters  aic  varied  and  the  problem  is 
solved  repeatedly  in  ‘he  forward  manner  until  a  best  fit  to  the  experiment  data  is 
obtained.  Figure  10  summarizes  the  optimum  parameters  for  the  model  to  give  the 
solution  shown  for  the  phase  velocity.  In  this  case,  the  best  fit  is  obtained  for  a  sediment 
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thickness  of  H  «  4.5  m  and  a  sediment  shear  velocity  vsi,  which  increases  linearly  with 
depth  z  as 

v$j  *  135  +  15  •  z  .  (2) 

The  most  important  parameters  are  layer  thickness  and  the  mean  shear  velocity  in  the 
layer. 


Frequency  -[HiJ 

Figure  10.  Bowidary  phase  velocity  as  a  function  of  frequency.  Theoretical  results 
calculated  with  H  =  45  m,  v^  =  135  +  15  xz,  pj/p0  *  2  and  P2/P0  -  3, 
v^  *  800  mis,  vcl  =  vc2  =  2000  mis  (from  Ref.  26). 

D.  Sediment  Shear  Properties 

Jensen  and  Schmidt29  give  simple  formulas  that  relate  the  speed  and  attenuation  cf 
interface  waves  to  those  of  body  shear  waves.  They  use  seismic  records  to  find  the  shear 
speed  and  attenuation  profiles  in  the  upper  sediment  layers. 

For  two  homogeneous  media  in  contact,  the  Scholte  wave  is  nondispersive.  Its 
propagation  speed,  Vsch,  and  attenuation,  Axfl,  are 
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-  0.9V, 

-ind 


(3) 


A^-l.lA,,  (4) 

where  V,  is  the  shear  velocity  and  A,  is  the  shear  attenuation  in  decibels  per  unit  length. 

These  formulas  are  accurate  to  within  a  few  percent  for  unconsolidated  sediments 
(clay,  silt,  and  sand)  with  low  shear  speeds.  Scholte  wave  properties  are  independent  of 
densities  and  compressional  wave  properties  in  the  two  media.30 

When  sediment  is  layered,  Scholte  waves  become  dispersive,  and  finite-ducted 
shear  modes  appear.  Figure  1 1  shows  dispersion  curves  for  a  sediment  layer  of  low  shear 
speed  p,  and  thickness  H  of  a  high-shear-speed  solid  acoustic  basement  or  substrate.31 
The  M’s  stand  for  the  mode  numbers.  The  lowest  mode,  M0,  is  an  interface  wave  that  in 


Figure  11.  Dispersion  curves  for  elastic  layer  on  a  semi-infinite  elastic  substratum 
(from  Ref.  29). 
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the  low-frequency  limit  kH  ->  0  propagates  as  a  Rayleigh  wave  on  the  substrate.  In  the 
high-frequency  limit  kH  -»  ••  the  lowest  mode  propagates  as  a  Rayleigh  wave  on  the 
surface  of  the  upper  layer.  The  surface  waves  all  have  a  cutoff  frequency  at  which  they 
propagate  with  the  shear  speed  of  the  substrate.  In  the  high-frequency  limit  these  modes 
propagate  with  the  shear  speed  of  the  upper  layer.  There  is  no  Stoneley  wave  because 
this  requires  a  change  in  shear  speed  of  <10%  across  a  solid-solid  boundary. 

Information  about  bottom  shear  properties  is  usually  obtained  from  measured 
dispersion  curves.  If  experiments  show  negligible  dispersion,  the  sediment  is 
homogeneous  to  a  depth  of  one  or  two  wavelengths,  which  is  the  penetration  depth  of  the 
Scholte  wave.  Equations  3  and  4  can  then  be  used  to  determine  the  shear  properties  of 
the  sediment. 

Previous  measurements  on  Scholte  waves  over  sand  or  silt,  between  3  and  35  Hz, 
have  shown  a  variation  of  attenuation  greater  than  two  or  more  orders  of  magnitude, 
from  about  0.02  to  2.3  dB/Xj.29  Here  X,  stands  for  the  shear  wavelength.  Very  low 
attenuation  values  go  with  low  frequencies  (3  to  5  Hz).  Two  possible  explanations  for 
this  large  variation  are:  (1)  shear  attenuation  decreases  with  depth  and  low-frequency 
Sholte  waves  penetrate  deeper,  and  (2)  attenuation  varies  with  a  power  of  frequency 
between  1  and  2.  The  Biot  model  predicts  such  a  dependence.32 

Two  types  of  measurements  can  be  used  to  determine  shear  speed  and  shear 
attenuation  profiles:  (1)  Stacked  time  signals  for  both  vertical  and  horizontal  particle 
velocities  are  shown  in  Figure  12  as  recorded  by  a  geophone  on  the  ocean  bottom.  The 
sources  were  explosive  charges  detonated  near  the  seafloor  in  20  m  of  water.  The 
dispersed  low-frequency  Scholte  waves  are  seen  to  arrive  with  group  velocities  between 
78  and  235  m/s.  (2)  By  using  a  multiple  filter  technique33  on  the  recorded  times  series  at 
a  fixed  range  of  1.7  km,  dispersion  diagrams  are  obtained  as  shown  in  Figure  13.  The 
contours  indicate  relative  energy  levels.  Note  that  there  is  more  information  in  the 
propagation  of  the  horizontal  radial  component  of  the  particle  velocities.  Figure  13 
shows  that  the  energy  arrives  in  three  discrete  modes  of  which  the  slowest  arrival  is  the 
Scholte  mode  (MQ)  with  energy  centered  at  about  2  Hz.  The  first  shear  mode  (Mj)  is 
strongly  excited  with  maximum  energy  around  2.8  Hz,  while  the  second  shear  mode  (M2) 
is  weak. 

The  measurements  are  then  applied  to  a  numerical  model  outlined  in  the  next 
section  in  a  trial-and-error  fashion  to  arrive  at  computed  dispersion  characteristics  that 
agree  with  the  measurements  of  Figures  12  and  13.  The  compressional  wave  properties 
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Figure  12.  Stacked  time  signals  for  vertical  (a)  and  horizontal  (h)  particle  velocities  as 
recorded  by  a  geophone  on  the  seafloor  (from  Ref.  29).  A  •  represents  one 
charge  unit. 
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Figure  13.  Dispersion  curves  obtained  by  applying  a  multiple  filter  technique  the 
experimental  records  at  range  R  =  1.7  km.  Energy  is  seen  to  be  propagat¬ 
ing  in  three  discrete  modes:  M0,  My,  and  M2  (from  Ref.  29). 

and  densities  are  fixed  a  priori  for  the  unconsolidated  sediment  since  they  are  known  to 
have  a  negligible  effect  on  the  propagation  of  interface  waves.  Figure  14  shows  the  final 
choices  of  the  shear  speed  profile  Cs  and  shear  attenuation  profile  Ps,  with  the  assumed 
compressional  wave  properties  Cc  and  Pr  and  density  p.  In  this  figure,  (3,  and  Pc  stand 
for  shear  and  compressional  wave  attenuation,  respectively. 
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Figure  14.  Model  environment  used  for  generating  synthetic  seismograms  (from  Ref. 
29). 


E.  Schmidt’s  Numerical  Method  j 

i 

Schmidt34,35  has  developed  a  fast  numerical  method  to  evaluate  viscoelastic 
physical  models  of  the  seabed  as  a  propagating  medium  for  seismic  waves.  He  analyzes 
the  excitation  and  propagation  of  interface  waves,  obtaining  many  of  the  basic  features 
observed  during  experiments  conducted  by  the  SACLANT  ASW  Research  Centre. 

One  of  the  purposes  of  investigating  interface  waves  is  to  determine  seabed 
properties  that  cannot  be  measured  directly.  A  distinct  transition  frequency  exists  below 
which  the  propagation  of  the  first  interface  wave  is  a  Rayleigh  wave  on  a  pure  acoustic 
basement  half-space.  Its  velocity  and  damping  yield  information  on  the  properties  of  the 
basement  The  transition  frequency  is  the  first  thickness-shear  frequency  of  the  sediment 
layer  as  an  infinite  plate  and  is  a  direct  measure  of  the  sediment  shear  velocity.  Above 
the  transition  frequency,  the  first  interface  mode  measures  the  properties  of  the  water- 
sediment  interface,  and  its  properties  determine  the  properties  of  the  sediment  layer.  The 
second  interface  mode  above  the  transition  frequency  is  related  to  the  sediment-basement 
interface  and  contains  coupled  information  on  both  media.  If  there  is  more  than  one 
sediment  layer,  several  interface  modes  will  appear;  their  interpretation  can  be  complex. 

The  data-fitted  theoretical  model  is  then  used  to  compute  dispersion  curves  and 
modal  energy  content  with  a  prescribed  frequency  dependence  of  shear  attenuation.  This 
dependence  turns  out  to  be  quadratic.  For  an  unconsolidated  sand-silt  bottom,  the  shear 
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velocity  is  found  to  increase  with  depth  from  85  m/s  at  the  seafloor  to  300  m/s  at  60  m 
depth.  The  shear  velocity  increases  with  depth  at  the  rate  3.58  m/s/m.  The  shear 
attenuation  at  3  Hz  decreases  with  depth  from  0.45  dB/X  at  the  seafloor  to  0.15  dB/X  at 
60  m  depth. 

Schmidt’s  technique35  is  similar  to  those  described  by  Kutschale36  and  DiNapoli 
and  Deavenport,37  who  take  a  Green’s  function  approach  based  on  the  Thomson-Haskell 
matrix  method.  However  Schmidt’s  method  is  computationally  mors  efficient. 

The  environment  is  approximated  by  a  series  of  horizontally  homogeneous  and 
isotropic  elastic  continua  with  Lame  constants  X  and  |i  and  density  p.  The  displacement 
and  stress  components  are  expressed  in  terms  of  three  scalar  potentials,  each  of  which 
must  satisfy  a  wave  equation.  After  the  Hankel  transform  is  applied,  the  boundary 
conditions  at  each  interface  yield  a  linear  system  of  equations  in  the  Hankel  transforms  of 
the  potentials.  A!  displacement  and  stress  components  can  then  be  expressed  as  Hankel 
transform  integrals  with  transforms  of  the  potentials  included  in  the  kernels. 

There  are  two  parts  to  the  numerical  evaluation.  First,  the  Hankel  transforms  of  the 
potentials  are  found  at  discrete  horizontal  wave  numbers.  The  equations  arc  then  solved 
conveniently  by  Gaussian  elimination.  Then,  the  desired  displacement  and  stress 
components  are  evaluated  at  a  given  depth  and  range  with  the  fast-field  integration 
technique  36,37 

Schmidt  applies  his  method  to  a  set  of  seven  test  cases  with  different  sediment 
materials  (silt  and  sand)  and  thicknesses  (5,  30,  50  m)  and  two  different  substrates 
(limestone  and  basalt).  He  also  studies  the  effect  of  a  linear  gradient  by  using  stepped 
layers.  The  water  depth  is  100  m  in  all  cases  with  a  source  amplitude  of  1  Pa  re  1  m  at 
50  m  depth.35 

He  uses  mode  transition  zones  from  the  theory  of  plate  vibrations,  where  they 
appear  as  thickness -shear  frequencies  for  an  infinite  elastic  plate  performing  free  shear 
vibrations  with  vanishing  vertical  displacements  at  the  edges.  For  a  silt  layer,  the  first 
thickness-shear  frequency  fTS  would  correspond  to  that  of  a  free  silt  plate  of  double  the 
actual  thickness 

/rs  =  CT/(4ds) ,  (5) 

where  CT  is  the  shear  velocity  and  ds  is  the  thickness  of  the  layer.  If  the  parameters  of 
silt  are  applied  to  Eq.  5,  then  we  get  fjs  =  1  Hz  (i.e.,  CT  =  200  m/s)  and  ds  -  SO  m. 
Thus 


fTS  =  200/200  =  1  Hz . 
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Because  the  transition  frequency  can  be  observed  experimentally,25  its  correlation 
with  the  thickness-shear  frequency  can  be  used  to  determine  the  shear  wave  velocity  in  a 
single  sediment  layer  overlying  a  rigid  half-space. 

In  general,  direct  measurements  of  shear  properties  are  difficult  to  perform. 
However,  both  shear  wave  speed  and  attenuation  profiles  can  be  inferred  from  measured 
dispersion  characteristics  of  bottom  interface  waves.  The  propagation  of  Scholte  waves 
depends  on  sediment  shear  properties.  Sediment  shear  wave  speeds  are  about  85  to  300 
m/s  with  shear  wave  attenuations  of  0.15  to  0.45  dB/X  at  3  Hz.  In  the  next  section,  I 
consider  the  remote  sensing  of  the  Biot  geophysical  properties  of  shear  modulus, 
permeability,  and  porosity,  which  allow  .  ,e  determination  of  the  geoacoustic  properties 
from  Biot  theory. 
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IV.  BIOT  GEOPHYSICAL  PARAMETERS 

A.  Introduction 

Yamamoto  and  his  associates  at  the  University  of  Miami  are  concentrating  on  a 
program  of  remote  sensing  and  inversion  of  bottom  seismic  data  to  determine  the  Biot 
geophysical  properties  of  marine  sediments.  They  are  particularly  interested  in  obtaining 
depth  profiles  of  porosity  and  permeability  from  synthetic  seismograms38  and  complex 
shear  modulus  profiles  in  the  sediments  as  a  function  of  depth  using  ocean  wave 
exciiation  of  the  bottom.6,39  Geoacoustic  properties  can  be  derived  from  these 
parameters.  This  unique  program  is  discussed  here. 

B.  Synthetic  Seismograms:  Porosity  and  Permeability 

Plane  wave  reflection  and  transmission  coefficients  play  an  important  role  in  the 
calculation  of  synthetic  seismograms  and  the  interpretation  of  acoustic  data. 

Claerbout40  showed  how  to  synthesize  a  layered  medium  from  its  acoustic 
transmission  response.  He  presented  a  direct  (noniterative)  inverse  method  of 
determining  an  acoustic  layered  medium  from  the  seismogram  caused  by  a  time-limited 
plane  wave  incident  from  the  lower  half-space  of  Claerbout’ s  model.  He  showed  that 
one  side  of  the  autocorrelation  of  the  seismogram  that  results  from  an  impulsive  source  at 
depth  is  the  same  as  the  seismogram  from  an  impulsive  source  on  the  surface.  This 
transforms  the  problem  to  the  acoustic  reflections  problem  solved  by  Kunetz.41  Both  the 
deep-source  time  function  and  the  layering  can  be  determined  from  a  surface 
seismogram. 

Goupillaud42  solved  the  inverse  problem  of  determining  the  medium  from  the 
seismogram.  Kunetz  showed  that  the  inversion  should  be  stable  and  will  be  so  if  the 
observed  reflection  seismogram  is  one  side  of  an  autocorrelation.  His  correlation 
function  is  the  autocorrelation  of  the  seismogram  observed  from  an  impulsive  source 
emerging  from  depth.  Thus  the  impulsive  transmission  problem  is  transformed  to 
Kunetz’s  reflection  problem  by  an  autocorrelation. 

Previously,  most  synthetic  seismograms  were  calculated  for  rigid  elastic  models  that 
ignored  the  effects  of  dispersion  and  attenuation.  In  a  more  realistic  model  of  marine 
sediments  in  which  the  elastic  material  is  replaced  by  viscoelastic  or  porous  viscoelastic 
material,  the  response  at  an  interface  between  two  different  layers  is  more  complicated. 
Reflection  and  transmission  coefficients  are  complex  numbers  for  a  viscoelastic  model. 
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Furthermore,  Deiesiewicz  and  Levy  have  shown  that  there  is  mode  conversion  between 
fast  and  slow  compressional  waves  and  generation  of  shear  waves  at  the  interfaces  of  a 

layered  porous  viscoelastic  medium.43 

Ganley44  outlined  a  method  that  allows  the  calculation  of  one-dimensional 
theoretical  seismograms  and  includes  the  effects  of  absorption  and  dispersion.  The 
parameters  that  are  specified  for  the  model  are  layer  thicknesses  and  density,  phase 
velocity,  and  the  specific  attenuation  factor  ( Q )  in  each  layer.  It  is  also  necessaiy  to  have 
equations  defining  the  absorption  and  dispersion  relations.  Futterman45  derived  a 
Kramers-Kronig  relation  that  allowed  the  calculation  of  dispersion  for  a  given  absorption 
function.  Ganley  used  one  of  Futterman’ s  absorption-dispersion  pairs  to  compute 
examples  at  the  surface  as  well  as  at  several  depths  for  a  simple  test  model. 

In  Ganley’s  model,  the  attenuation  is  accounted  for  mathematically  by  allowing  the 
elastic  modulus  to  be  a  complex  function  of  frequency  (viscoelastic  model).  This  results 
in  a  complex  velocity  and  wave  number,  the  reflection  and  transmission  coefficients  also 
become  complex  functions  of  frequency.  The  attenuation  model  is  applicable  to  plane 
waves  in  a  flat,  layered  model. 

Turgut  and  Yamamoto38  use  a  more  realistic  Biot  model  of  marine  sediments  than 
Ganley  in  which  the  viscoelastic  material  is  replaced  by  a  porous  viscoelastic  material  for 
which  the  response  at  an  interface  between  two  different  layers  is  more  complicated.  In 
this  latter  case,  Biot  theory  predicts  mode  conversion46,47  between  fast  and  slow 
compressional  waves  at  interfaces.  Their  model  includes  this  mode  conversion  as  an 
energy  loss  mechanism.  This  type  of  energy  loss  is  even  more  important  when  multiple 
reflections  occur  within  layers.  They  present  a  possible  method  of  predicting  porosity 
and  permeability  of  marine  sediments  by  using  information  on  the  frequency  dependence 
of  the  phase  velocity  and  attenuation  of  the  fast  compressional  waves.  Biot  theory 
enables  the  determination  of  porosity  and  permeability  of  marine  sediments  once  the 
frequency-dependent  velocity  and  anelasticity  iQ'1)  are  measured. 

1 .  Determination  of  Porosity  and  Permeability 

The  propagation  of  acoustic  waves  through  marine  sediments  can  be  considered  as  a 
fluid-solid  (inertial  and  viscous)  coupling  problem  governed  by  hydrodynamic  properties 
of  sediment,  porosity,  and  permeability.  The  inertial  coupling  is  governed  by  porosity; 
the  viscous  coupling  by  permeability.  Jackson  and  Anderson  point  out  that  a  viscous 
coupling  is  a  relaxation  process.48  Therefore,  the  velocity  V  and  anelasticity  Q~l  for  the 
fast  compressional  wave  are  frequency-dependent  Turgut  and  Yamamoto  calculate  the 
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compressional  wave  velocity  and  anelasticity  from  the  Biot-Stoll  equation  using  the 
permeability  ks,  porosity  |3,  shear  modulus  G,  and  frame  loss  coefficient  5  of  coarse  sand. 
Figure  15  shows  that  the  inelasticity  becomes  maximum  at  about  32  Hz  for  coarse  sand. 
This  peak  frequency  is  the  relaxation  frequency  of  the  sand.  The  relaxation  frequency  of 
a  saturated  sediment  is  given  by  Yamamoto44  as 

f  r  =  0. 1  -j~- ,  (6) 

"j 

where  v  is  the  kinematic  viscosity  of  pore  fluid. 


Figure  15.  Propagation  velocity  VR  and  anelasticity  Q'1  values  for  coarse  sand  (from 
Ref.  38). 


The  weak  frame  (Kr,  Kj  »  Kb  p)  approximations  in  Plona  and  Johnson49  are  valid 
for  most  marine  sediments.  Kr  is  the  bulk  modulus  of  the  grains,  Kj  is  the  bulk  modulus 
of  the  fluid,  Kb  is  the  bulk  modulus  of  the  frame,  and  p  is  the  shear  modulus  of  the 
frame.  Thus,  we  can  use  the  following  approximate  relations  in  their  respective 
frequency  bands. 
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a.  Low-frequency  range 


VmVmmlMi°Lr 


and 


Q-l-  271--  —f+  5'  , 

PP/  V 


(7) 


where  V  is  the  compressional  wave  velocity  and  subscript  o  stands  for  low-frequency 
limit. 


8'=  50' 


pv, 


where  5  is  the  skeletal  frame  loss,  and 


2  * 


2  n 


G'  =  G(2  +  -~~). 

1  -n 

Here  n  is  the  Poisson  ratio  and  G  is  the  shear  modulus. 


b.  Near  relaxation-frequency  range 


G-1  =  e^,=yP-^£-|^  +  »'.  (8) 

3  P  P/ 

If  we  measure  V  and  Q1  of  a  sediment  as  a  function  of  frequency,  we  can  determine 
the  porosity  and  permeability  from  equations  7  and  8.  In  addition,  the  frame  loss  5  can 
be  determined  if  the  shear  modulus  G  is  known. 

Figure  16  from  Reference  38  shows  synthetic  seismograms  calculated  for  a  Ricker 
w  .velet  propagating  through  homogeneous  layers  of  (a)  silt,  (b)  medium  sand,  and  (c) 
arse  sand.  A  Ricker  wavelet  has  the  spectrum 

N  (to)  =  co2  exp(-co2 / co2) ,  (9) 

where  co  is  the  angular  frequency  and  co0  is  the  peak  angular  frequency. 
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C.  Sediment  Shear  Modulus 

Yamamoto  and  Torii6  have  developed  a  theory  that  enables  seabed  shear  modulus 
profiles  to  be  extracted  by  inverting  measurements  of  seabed  motion  and  pressures 
induced  by  water  waves  at  one  point. 

As  surface  water  waves  propagate  through  intermediate  to  shallow  water,  they 
induce  motions  of  the  seabed  sediments  through  variations  in  pressure  on  the  seafloor. 
According  to  models  developed  by  Yamamcto  and  others,  the  layered  seabed  sediments 
can  be  seen  to  behave  in  a  massless,  incompressible,  elastic  manner  in  response  to 
pressures  induced  by  water  waves.50-52  Under  these  circumstances  one  can  predict 
realistically  the  seabed  response  to  passing  water  waves  if  the  sediment  shear  modulus  is 
known  at  every  depth.53  Conversely,  through  the  use  of  geophysical  inverse  methods, 
one  can  extract  the  profile  of  the  shear  modulus  as  a  function  of  depth  from 
measurements  of  seabed  motion  and  pressure  induced  by  water  waves. 

Linear  inverse  theory  is  used  combined  with  Yamamoto’s  theory  of  wave-seabed 
interaction  (see  Appendix  A).  The  uniqueness  and  consistency  of  the  inversion  was 
established  by  numerical  tests  using  synthedc  data  without  noise.  As  a  practical  test  on 
noisy  data  it  was  shown  that  good  agreement  was  obtained  between  the  inverted  shear 
modulus  profile  and  direct  measurements  at  the  Mississippi  River  Delta.  Theoretical 
assumptions  of  small  amplitudes  and  elastic  sediment  behavior  were  validated  by 
examining  time  series  of  seabed  displacement,  dynamic  phase  spectra,  and  profiles  of 
seabed  shear  strain. 

Based  on  the  success  of  these  experiments,  a  Bottom  Shear  Modulus  Profiler 
(BSMP)  was  developed.54  Comparisons  of  BSMP  inversion  results,  in  situ  sediment 
strength  tests,  and  laboratory  analysis  of  core  samples  show  that  this  technique  can 
accurately  measure  the  magnitude  and  depth  structure  of  the  sediment  shear  modulus. 
The  shear  modulus  of  a  marine  sediment  is  one  of  the  more  important  parameters  for 
many  theories  modeling  seismic,  acoustic,  and  surface  gravity  wave  propagation  in  the 
ocean.  Existing,  accepted  methods  of  measuring  a  profile  of  the  sediment  shear  modulus, 
such  as  making  standard  penetration  boreholes  or  measuring  rhe  travel  time  of  cross-hole 
shear  waves,  are  cumbersome  and  slow.  Such  methods  also  suffer  from  the  unavoidable 
problem  that  the  act  of  boring,  drilling,  or  removing  the  sediment  disturbs  the  sediment 
framework  and  grain  structure.  This  can  seriously  affect  the  sediment  shear  modulus, 
especially  for  cohesive  muds  and  clays.  In  contrast,  this  new  BSMP  system  is  fast  and 
unobtrusive  to  the  seabed.  It  seems  possible  that  the  BSMP  could  also  become  a 
standard  tool  for  engineering  surveys  of  marine  sediments. 
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Additional  experiments  were  conducted  off  the  New  Jersey  coast  in  August 
1986.54,55  Of  concern  was  the  fact  that  the  accurate  measurement  of  the  small  seabed 
motions  occurring  in  a  typical  marine  environment  is  not  a  trivial  problem.  The  typical 
seabed  displacement  amplitude  caused  by  a  moderate  sea  on  a  sandy  bottom  is  of  the 
order  of  a  millimeter.  A  further  problem  is  the  extremely  low  frequency  of  surface 
gravity  water  waves,  which  are  typically  in  the  range  of  0.06  to  0.3  Hz  (periods  of  17  to  3 
seconds).  The  motion  sensors  need  to  be  firmly  coupled  to  the  seabed,  allowing  no 
slippage  or  relative  motion  between  the  sediments  and  the  instruments.  This  seabed-to- 
instrument  coupling  is  a  difficult  quantity  to  assess  and  may  be  a  significant  problem  on 
softer,  more  slippery  muds  and  clays.  Also,  the  instrumentation  package  as  a  whole 
needs  to  be  waterproof  and  pressure-resistant  to  depths  of  200  m  or  more.  Figure  17 
shows  a  detailed  diagram  of  the  major  components  of  the  BSMP  instrumentation 
package.  The  main  body  of  the  instrument  housing,  the  seismometer  supports,  and  the 
bottom  base  plate  are  made  of  aluminum  alloy.  The  hemispherical  glass  cap  is  designed 
to  withstand  a  depth  of  600  m  and  to  have  a  low  profile  in  order  to  minimize  horizontal 
hydrodynamic  forces  on  the  unit.  In  use,  a  slight  vacuum  is  applied  inside  the  glass 
hemisphere  to  hold  the  cap  in  place  and  reduce  leakage  by  compressing  the  "O"  rings. 
Fully  assembled,  the  entire  BSMP  unit  has  a  mass  of  about  45  kg. 

The  heart  of  the  BSMP  sensor  is  formed  of  three  seismometers  mounted  in  one 
vertical  and  two  perpendicular  horizontal  directions.  Each  seismometer  measures 
motions  along  its  principal  axis  only,  and  the  three  directions  are  labeled  vertical, 
horizontal,  and  transverse.  A  compass  and  pendulum  tiltmeter  are  included  in  each 
BSMP  unit  to  determine  the  orientation  of  the  seismometers.  The  sensitivity  of  the 
seismometers  is  45  volts/m/s  with  a  resolution  threshold  of  3  x  10'8  m/s. 

Figure  18  shows  the  sensitivity  and  frequency-corrected  BSMP  time-series  for  the 
four  measured  quantities  required  for  the  inversion  process.  Figure  19  shows  the  velocity 
and  pressure  power  spectra  indicating  the  ocean  wave  (swell)  periods  contributing  energy 
for  this  technique.  In  Figure  20,  the  important  parameters,  admittance  and  coherence,  are 
plotted  for  one  experiment  Figure  21  provides  a  comparison  of  the  shear  modulus 
inversion,  reference  data  from  a  sediment  shear  strength  borehole,  and  laboratory 
analysis  of  a  vibracore  sample.  The  agreement  is  remarkable. 
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Figure  17.  The  waterproof  housing  and  various  instruments  that  comprise  the  Bottom 
Shear  Modulus  Profiler  (from  Ref.  39). 
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RADIAL  SEISMOMETER  TRANSVERSE  SEISMOMETER 


Figure  19.  Frequency-corrected  and  time-averaged  seismometer  velocity  and  pressure 
power  spectra  from  experiment  121  _5  (from  Rtf.  39). 
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Figure  20.  Vertical  and  radial  admittance  data  with  associated  coherences  versus 
wave  period  from  experiment  121  _5  (from  Ref.  39). 
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Figure  21 .  Comparison  of  inverted  shear  modulus  profile  from  experiment  121  5  with 
reference  data  from  a  sediment  shear  strength  borehole  and  laboratory 
analysis  of  a  vibracore  sample.  The  inversion  used  only  vertical  admit¬ 
tance  data,  and  an  eigenvalue  expansion  limit  of  5  (from  Rtf.  39). 
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V.  CONCLUSIONS 

The  methodology  of  geoacoustic  and  geophysical  inversion  techniques  has  now 
been  developed  and  proven  to  the  point  that  acoustical  experiments  and  surveys  in 
shallow  water  should  use  it  on  a  regular  basis.  When  the  inversion  met  lodology  is 
combined  with  remote  sensing  instrumentation,  we  have  a  fast  and  efficient  way  to  map 
the  geoacoustic  properties  of  the  seabed  as  a  function  of  depth. 

For  shallow  water  areas,  the  perturbative  inverse  technique  based  on  measurements 
of  the  depth-dependent  Green’s  function  is  adequate  to  construct  the  sound  speed  and 
attenuation  profile  of  the  compressional  waves.  Synthetic  seismograms  using  a  Biot 
geoacoustic  model  can  be  used  to  obtain  interface  losses  resulting  from  compressional 
waves  as  well  as  porosity  and  permeability. 

The  study  and  measurement  of  seismic  interface  waves  in  the  dispersive  range  of 
0. 1  to  0.025  s  provide  information  on  the  shear  wave  properties  of  the  bottom  In 
addition,  it  has  been  shown  that  ocean  gravity  swell  waves  of  period  5  to  10  s  provide 
data  for  the  construction  of  accurate  shear  modulus  profiles  as  a  function  of  depth. 
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Appendix  A 


Summary  of  Inverse  Theory  for  Bottom  Shear  Modulus  Profiles 


* 


This  present  summary  has  been  adapted  from  Yamamoto  and  Torii  A_1  with  the  authors’ 
kind  permission. 
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I.  UNEAR  INVERSE  THEORY 

The  linear  inverse  problem  most  often  encountered  in  geophysics  can  be  formulated 
as 


*i  =J i  Ci(y  )m  (y  )dy  (i  =  1,2,3,.., N) 

or,  in  the  computationally  simpler  discrete  form, 

M 

=  £  Cijntj  (i  =  1,2, 3...., AO . 

>=i 

In  matrix  notation 


(A-l) 


(A-2) 


e  =  Cm  ,  (A-3) 

where  et  is  one  of  N  observations  (DATA),  C  is  a  kernel  derived  from  theory,  and  m  is 
the  unknown  function  (MODEL)  defined  on  the  I  interval.  A,2“A4 

The  singular  value  decomposition  of  matrix  C  is  l  "d  to  solve  equation  (A-3) 
following  Lanczos^-5 

C  =  UAVt  .  (A-4) 

[/,  :  the  columns  of  U  are  the  N-dimensional  eigenvectors  of  a  coupled  eigenvector- 
eigenvalue  problem  that  spans  the  data  space.  V,  :  the  columns  of  V  are  the  M- 
dimensional  eigenvectors  of  a  coupled  eigenvector-eigenvalue  problem  which  spans  the 
solution  (MODEL)  space.  A  is  the  diagonal  matrix  with  the  corresponding  eigenvectors 
X;,  which  are  sorted  into  descending  order.  Zero- valued  eigenvectors  are  not  included. 
This  leads  to 

CCTUi  =  XfUi  (i  =  l,2 . N) , 

C  TCVi  =  Xfy  (i  =  1 ,2 . M ) ,  (A-5) 

and  the  C7,  and  V,-  are  orthonormal  such  that 

VjUj-hij,  VjVj-Sij.  (A-6) 

In  terms  of  U,  V,  and  A,  the  solution  to  equation  (A-3)  is 
m  =C~le 


=  VA ~lUTe 
N  U:  '  e 

=  Z  -\r  V*'  • 

1=1  A* 
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For  the  data  with  noise,  the  following  is  considered  for  the  construction  of  the 
solution  using  the  truncated  singular  value  expansion, 

A  L  Ur  e 

**  *  E  -T7-  Vi  (L<N) .  (A- 8) 

i=i  M 

There  are  good  reasons  for  using  the  truncated  expansion.  With  the  assumption  of 
Gaussian  white  noise,  one  can  show  that  the  uncertainty  of  each  coefficient  multiplying 
V,  in  equation  (A-7)  is  proportional  to  X,-1/2.  Thus,  the  expansion  equation  (A-7)  is  in 
terms  of  vectors  whose  coefficients  increase  in  uncertainty.  It  is  also  usually  the  case  that 
the  vectors  Vf  become  more  oscillatory  as  i  increases  (exceptions  to  this  rule  are  very 
rare).  Thus,  the  truncated  expansion  equation  (A-8)  simultaneously  removes  those  parts 
of  the  solution  that  are  sensitive  to  noise  and  gives  a  smoothed  estimate  of  the  solution. 
This  point  is  demonstrated  here  using  synthetic  data  without  noise  and  measured  data 
with  noise. 

The  resolution  matrix  has  been  adopted  to  determine  the  quality  of  the  above 
solution  using  the  truncated  singular  value  expansion 

c-c-U-v,.  (A-9) 

In  the  ideal  case  in  which  C_1C  results  in  the  identity  matrix,  the  resolution  is  perfect. 
When  C~lC  has  near  unity  values  on  the  diagonal  and  small  values  elsewhere,  the 
computed  solution  represents  a  smoothed  version  of  the  true  solution. 

H.  NONLINEAR  INVERSE  PROBLEM 

All  of  the  above  analyses  apply  to  the  linear  inverse  problem.  For  the  nonlinear 
problem  such  as  the  one  under  consideration,  an  iterative  linearization  scheme  has  been 
used  by  which  successive  approximations  to  the  model  can  be  computed .  Starting  with 
the  nonlinear  relation  between  data  e  and  model  m ,  and  recalling  equation  (A-3), 

e  =/(m) .  (A- 10) 

In  order  to  minimize  4>  =  |  e  —f  (m  0)  | 2,  it  is  necessary  to  assume  an  initial  value  m  0  and 
then  linearize.  Expanding  equation  (A- 10)  in  a  Taylor  series  form  results  in 

e  =/(m)=/(mo)+  (m  -m  0)  (A- 11) 

om  m=m0 

or 
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<  ~/(™o)  =  (-^4^)  (m-m o) 

dm  '  mmn0 

5(/(m*H-/W)  (A-12) 

Am  '»»■«»• 

The  derivatives  df  /dm  are  approximated  by  finite  differences  in  the  last  equation.  As 
Yamamoto  shows,  the  finite  difference  approximations  to  the  derivatives  are  quite  good 
for  our  problem.  His  model  and  parameters  allow  him  to  arrive  at  fire  following  matrix 
equation: 

Ac  =  C  Am  .  (A- 13) 

The  linear  system  expressed  by  equation  (A- 13)  is  solved  for  Am  using  equation  (A-8) 
and  then  repeating  the  process  after  replacing  m0  by  m0  +  Am .  This  process  is  repeated 
until  Am  becomes  acceptably  small.  In  practice  the  method  does  not  always  converge. 
A  serious  shortcoming  of  this  linearization-iteration  approach  is  that  even  when  a 
solution  is  found  it  is  not  possible  to  say  without  testing  whether  another  equally 
acceptable  (but  well  separated  in  model  space)  model  exists. 

When  problems  are  treated  as  nonlinear  it  is  found  that  only  the  second,  third,  or 
fourth  terms  are  retained  for  the  first  few  iterations,  but  when  approaching  the  final 
solution  the  solution  retains  more  terms. 
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